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Abstract 


López-Acosta, N. P., £ González-Acosta, J. L. (September- 
October, 2015). Study of Water Flow in Dams using 
Successive Over-relaxation. Water Technology and Sciences (in 
Spanish), 6(5), 43-58. 


Free surface problems represent boundary value problems in 
which a portion of the boundary, the free surface, is unknown 
and must be determined as part of the solution. Classical 
rigorous and approximate procedures used to draw this 
line are limited to homogeneous and isotropic media with 
specific geometries. Currently, this can be determined using 
numerical methods such as the finite element method (FEM). 
Nevertheless, FEM requires the storing and handling of a 
large number of matrices to solve linear equation systems, 
increasing calculation time. The present article proposes an 
alternative to analyze free surface problems based on the 
numerical solution of finite difference equations using the 
successive over-relaxation method (SOR). Two techniques 
are implemented with the SOR method— Baiocchi's Solution 
and the Extended Pressure Method with the iterative Gauss- 
Seidel process. First, the theoretical basis for these methods 
are provided. Then, their applicability is described according 
to an analysis of unconfined flow in a homogeneous and a 
heterogeneous dam. As part of the results, the upper flow 
lines obtained with each technique and the flow networks 
calculated with the SOR method are presented and the use 
of finite difference equations to determine the hydraulic 
gradient and flow rate through the flow domain is described. 
Lastly, conclusions and recommendations for applying and 
optimizing the use of these techniques are provided. 


Keywords: Free surface, unconfined flow, Successive Over 
Relaxation (SOR), Gauss-Seidel's iterative process, finite 
differences, upper flow line, Extended Pressure method, 
Baiocchi's method, flow networks in dams. 


Resumen 


López-Acosta, N. P, € González-Acosta, J. L. (septiembre-octubre, 
2015). Estudio del flujo de agua en presas con sobre relajaciones 
sucesivas. Tecnología y Ciencias del Agua, 6(5), 43-58. 


Los denominados problemas de superficie libre son problemas de 
valores de frontera en los que una porción de la frontera, la superficie 
libre, no se conoce y debe determinarse como parte de la solución. 
Los procedimientos rigurosos y aproximados clásicos para el trazo de 
esta línea se limitan a medios homogéneos e isótropos con geometrías 
particulares. Actualmente es posible emplear métodos numéricos 
como los elementos finitos (MEF) para su determinación. Sin 
embargo, el MEF requiere la resolución de sistemas de ecuaciones 
lineales que involucran el almacenamiento y manejo de un gran 
número de matrices que incrementan el tiempo de cálculo. En este 
artículo se propone una alternativa para analizar problemas de 
superficie libre que se basa en la resolución numérica de ecuaciones de 
diferencias finitas mediante el método de sobre relajaciones sucesivas 
(SOR, Successive Overrelaxation). Se expone la implementación de 
dos técnicas basadas en el método SOR: la Solución de Baiocchi y 
el Método de la Presión Extendida, aplicando el proceso iterativo 
de Gauss-Seidel. Inicialmente se proporcionan las bases teóricas de 
estos métodos. De manera posterior, su aplicabilidad queda expuesta 
con el análisis del flujo no confinado en una presa homogénea y en 
otra heterogénea. Como parte de los resultados, se presentan las 
líneas de flujo superior obtenidas con cada técnica y las redes de flujo 
calculadas con el método SOR; asimismo, se explica cómo determinar 
el gradiente hidráulico y el gasto de infiltración a través del dominio 
de flujo mediante ecuaciones de diferencias finitas. Por último, se 
proporcionan comentarios concluyentes y recomendaciones de cómo 
aplicar y optimar el empleo de estas técnicas. 


Palabras clave: superficie libre, flujo no confinado, Sobre 
Relajaciones Sucesivas (SOR), proceso iterativo de Gauss-Seidel, 
diferencias finitas, línea superior de flujo, Método de la Presión 
Extendida, método de Baiocchi, redes de flujo en presas. 
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Introduction 


The classical relaxation method is an iterative 
process in which solutions for water flow 
through porous media can be obtained by 
simply knowing the domain geometry and 
hydraulic boundary conditions (Allen, 1954; 
Finnemore éz Perry, 1968). This method can solve 
Laplace's equation for a point (node) relative to 
its surrounding points using an algebraic finite 
difference equation (Southwell, 1940). For this 
procedure, a square mesh with dimensions of 
A, = A, is drawn in the flow zone if the medium 
is homogenous and isotropic, and a rectangular 
mesh with dimensions of A, + A, is drawn if the 
medium is anisotropic. The intersections of the 
squares or rectangles constitute the nodes of the 
mesh. For these nodes, approximate values of 
the hydraulic head or potential / (points where 
h requires to be calculated) must be assigned 
while respecting the known values of h in the 
flow boundaries. These values usually cor- 
respond to the upper and lower water levels 
or the upstream water level (UWL) and down- 
stream water level (DWL) of the problem at 
hand (known as the Dirichlet boundary conditions; 
Cheng 4 Cheng, 2005). The values assigned in 
the nodes are arbitrary and can be zero or the 
result of a reasonable estimation. However, 
although there are several techniques that can 
be used to ensure that the value of the potential 
imposed on the nodes where h is not known 
is as accurate as possible (Young, 1950), it is 
important to verify the precision of the assigned 
data manually by calculating the residue in each 
node. For example, the difference between the 
hydraulic potential of the four surrounding 
nodes is calculated with regard to the central 
or interior node, and so on. Therefore, the 
relaxation procedure involves the systematic 
refinement of this residue throughout the grid 
until the residue in all mesh nodes of interest 
is zero or practically zero. This value indicates 
compliance with Laplace's equation when ex- 
pressed as a finite difference and results in a 
solution for a certain water flow problem (Allen, 
1954; Finnemore éz Perry, 1968). The classical re- 


laxation method has several disadvantages. For 
example, this method becomes long and tedious 
when trying to obtain approximate solutions if 
the method has not been programmed using a 
computer. Besides, this method does not pro- 
vide general solutions and must be applied to a 
particular case. However, the classical relaxation 
method can be adapted for diverse conditions. 
Due to the versatility that is demonstrated for 
this method, several improvements have been 
made, including optimization of the calculation 
time, taking into account the variations of the 
materials in the flow region, and solving more 
complex problems, such as free surface prob- 
lems (unconfined flow). 

The Successive Overrelaxation (SOR) method 
(Young, 1950) is a modification to the classical 
relaxation method. This method is powerful 
and can be used to obtain approximate nu- 
merical solutions for multiple equations with 
unknown analytical solutions (Young, 1950). 
An important advantage of the SOR method 
is that it automatically solves the variations in 
the hydraulic potential for all mesh nodes that 
represents the flow region. The SOR method 
uses the Gauss-Seidel iterative process, which 
utilizes cyclic finite difference equations to 
solve the problem of interest (Wang $ Ander- 
son, 1982). Applying this method, the solution 
process is automatically repeated for each node 
until the desired tolerance or minimum error 
is obtained. For unconfined flow problems in 
which the upper flow line or free surface is un- 
known a priori, methods have been devised to 
modify the ruling equations of the nodes. These 
methods can be used for solving free surface 
problems with the basics of SOR method, such 
as the Baiocchi's solution (Baiocchi, 1971; Bruch, 
1980) and the Extended Pressure method (Brezis, 
Kinderlehrer, €: Stampacchia, 1978; Bardet éz 
Tobita, 2002). The first solution is useful for 
homogenous media, and the second solution 
is useful for homogenous and heterogeneous 
media. However, these methods have not been 
applied to practical cases. 

Generally, free surface problems represent 
a Challenge in many areas of fluid mechanics 
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because they involve boundary value problems 
in which a portion of the boundary is unknown 
and must be determined as part of the solution 
(Cryer, 1970). The presence of the free surface or 
water table makes the analysis methods more 
difficult. Dupuit's parabola (Dupuit, 1863) and 
Kozeny's parabola (Kozeny, 1931) are rigorous 
solutions for drawing the upper flow line and 
are only applicable for homogenous and isotro- 
pic media with vertical walls (Dupuit) or with 
filters (Kozeny). Other approximated methods, 
such as the tangent (Schaffernak, 1917; Iterson, 
1916, 1917) and sine methods (Casagrande, 
1932), allow mainly calculate the discharge 
point of the upper flow line. Currently, numeri- 
cal methods, such as the finite element methods 
(FEM), can be used to determine the discharge 
point of the upper flow line. However, the FEM 
requires solving a system of linear equations 
and involves storing and handling matrices that 
increase the calculation time. 

This paper aimed to implement the SOR 
method as an alternative for analyzing free 
surface or unconfined flow problems in ho- 
mogenous and heterogeneous media. First, the 
theoretical background of water flow in soils is 
provided, especially the solution to Laplace”s 
equation using finite differences. Next, the basic 
equations of the SOR method are obtained. Sub- 
sequently, this technique is enabled to evaluate 
unconfined flow problems using the Baiocchi's 
and Extended Pressure methods. The applica- 
bility of these techniques is demonstrated by 
analyzing the unconfined flow of water through 
a homogenous earth dam and a dam composed 
of different materials. We explain how to deter- 
mine the upper flow line and use its position to 
solve the flow problem by considering confined 
flow and modifying the boundary conditions 
and the fundamental equations to determine the 
equipotential and flow lines. Thus, the upper 
flow line is obtained for each technique and the 
flow nets are calculated with the SOR method. 
In addition, we explain how to calculate the 
hydraulic gradient and the seepage flux or flow 
rate through the flow domain using finite differ- 
ence equations. Finally, we present concluding 


remarks and recommendations regarding the 
application of the SOR method. 


Theoretical Basis of Water Flow in Soils 
Laplace's Equation 

The laminar flow through porous media obeys 
Darcy's law and is governed by the following 
equation: 


v= k—= (1) 


where: 
v = discharge velocity. 


total pressure or total hydraulic head. 


h 
k = permeability of the soil. 

Il = distance traveled by a water particle. 
In vector format (Harr, 1962): 


v=-kgrad(h) (2) 


Furthermore, the total pressure h is defined 
by the Bernoulli's law as follows (neglecting the 
velocity head): 


h=[(p/p3)+y] (8) 


where: 

p = hydrostatic pressure (or pore pressure). 

p = water density. 

g = acceleration due to gravity. 

y = position pressure with respect to a level of 
reference. 


Because the amounts of incoming and outgo- 
ing flow in the normal direction n to the faces 
of a cubic element with dimensions of dx, dy 
and dz is the same (due to the principle of flow 
continuity), the following equation holds: 


a das ds 


9x dy 


A ai (4) 
OZ 
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By factoring and removing terms, the conti- 
nuity equation is obtained in three dimensions 
as follows: 

dv, 00, ¿vu 
— 424 


—Ñ=0 5 
Ox  0y  0z 6) 


To solve flow problems, it is convenient to 
use a velocity potential function q rather than the 
total pressure h. The q is described as follows: 


d(x,y,z) =-kh (6) 


Therefore, the discharge velocity in each of 
the main directions becomes: 


Thus, Laplace's equation is obtained by substi- 
tuting equation (7) in equation (5): 


2 2 2 
0 0, (8) 
dx” 0y”  0z 


If the medium is homogenous and isotropic 
and flow only occurs in two dimensions (X and 
Y, two-dimensional flow), then: 


Se, Po (9) 


Approximate Solution of Laplace's Equation 
by Means of Finite Differences 


The basic concept of the finite difference method 
is based on replacing the continuous partial de- 
rivatives of the equation that governs water flow 
with the variation ratio of the x and y variables 
in small finite increments (Mahmud, 1996). In 
this way, the following equations are obtained 
when applying Taylor series expansion to the 
potential function q (Mitchell € Griffiths, 1980): 
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0 1%,» 
Piaj=0,,+ a + e (10) 
and: 
0% 19% .2 
dr; =0,; ro (11) 
where: 


Ax, = distance in the X, direction between two 
nodes. 


Ax, = distance in the X, direction between two 
nodes. 
row in which the node of interest is lo- 
cated. 


a. 
Il 


j  = column in which the node of interest is 


located. 


Figure 1 shows the representation of a node 
that is affected by four surrounding nodes and 
can be used to observe the previously men- 
tioned variables. 

Expression (10) is an approximation to the 
solution of Laplace's equation regarding the 


central and the frontal nodes (dp ) Equation 


i+l, 
(11) represents an approximation to the solution 
of Laplace's equation regarding the central and 
rear nodes (d,,)- The equations for the upper 
(d 1) and lower (O, ;.1) nodes are analogues to 
expressions (10) and (11). From equations (10) 
and (11), and by applying a series of mathemati- 
cal manipulations, the solution of the potential 
function q, is given by: 


ij 
7 
Ay 
Vi db. do 
-1/ 1, i+1 
id e e 
Ay, 
jj 
5d 
Ar, Axo 


Figure 1. Arrangement of nodes in a flow region. 
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Ax? 
Deja + Pr + yr ia E di: ,) 
b,,= . (12) 
ij LE 
q 
Ay 
2 
I£A ==, then: 
CO FRA O +0. 
d,, ES Deje Pon (9... di) (13) 


2(1+12) 


Equation (13) approximates the solution 
of the Laplace's equation by the classical re- 
laxation method in terms of finite differences. 
From this approximation, the solutions of the 
nodes located inside the flow region of interest 
are obtained. Likewise, the expressions for the 
nodes in other locations (impervious boundar- 
ies, corners, etc.) are determined. 


Basis of the Basic Equations of the SOR 
Method 


From equation (13), it is possible to apply a 
series of algebraic operations to obtain the ba- 
sic equations of the successive overrelaxation 
method by assuming the following (Young, 
1950): 


01), aa) 

where: 

Cc  =residue after two Gauss-Seidel iterations 
in a single node. 

m  =iteration number. 

d, = value of the potential function in the node 
1,| in iteration m. 

e = value of the potential function in the same 
node ¿,¡ in iteration m+1. 


By solving 4/7 from (14), the following equa- 
tion is obtained (Young, 1950): 


9 =0)+ac (15) 


where ais a relaxation factor thatis added to the 
equation. For fast convergence, it is recommend- 
ed that the relaxation factor is within the inter- 
val of 1 < a.<2 (Young, 1950). In fact, the value 
of a, is used to name this method as follows: 
over-relaxation if 1 <a <2, and under-relaxation if 
0 <a <1 (for a > 2 the process diverges). Ad- 
ditionally, if equation (13) is substituted in equa- 
tions (14) and (15), then (Wang éx Anderson, 
1982): 


e =(-a)jp, 


jar 0 0 (0%, +003)) 
a (16) 
2(1+27) 
Expression (16) represents the solution 

to Laplace's equation for a central node with 
respect to the adjacent nodes using the SOR 
method. The expressions described with this 
method for evaluating the potential function y 
(Wang € Anderson, 1982) are applicable to other 
flow conditions that are discussed in this paper. 
Thereby, for a node located at the boundary of 
an impervious horizontal surface, the solution 
to Laplace's equation with the SOR method is: 


1 1 
7 =(1-aJo"; 
Ql 


ee) ate] 0 


For a node located in a vertical impervious 
boundary, the solution is: 


a 


tem) a] a 


For the nodes located in the vertices or exte- 
rior corners, the solution is: 


bs a (a EN 0; 


DRA) + (20, +07) 


3(1+22) Je 
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For the nodes located in the transition zone 
of two materials with different permeability's, 
the solution is: 


e = (1-9; 


0002 Jarre) een, 


2 1+y 


+ (20) 


For a node located in the transition zone 
of two materials and at an impervious lower 
boundary, the solution is: 


e = 0-0, 


m+1 m+1 m+1 m 
q v(9r7, +), +; (21) 
2 1+y 


where y = k,/k,, k, is the permeability of the first 
material and k, is the permeability of the second 
material. 


Implementing the SOR Method for 
Problems of Unconfined Water Flow 


If the zones of the flow domain where the water 
pressure is zero are determined for free surface 
or unconfined flow problems, it is assumed that 
the upper flow line is defined at these points 
(Cryer, 1970). The Baiocchi's solution (Baiocchi, 
1971; Bruch, 1980) and the Extended Pressure 
method (Brezis et al., 1978; Bardet €: Tobita, 2002) 
are two variants of the SOR method that can be 
used to determine the position of the upper flow 
line in homogenous media or media composed 
of materials with different permeability values, 
respectively (as explained below). 


Solution to the Problems in Homogenous 
Media Using the Baiocchi's Method 
(Baiocchi, 1971; Bruch, 1980) 


Instead of using the velocity potential function 
q, the Baiocchi's method mathematically solves 
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Laplace's equation (in the flow zone and in its 
boundaries) for unconfined flow by using a 
new variable 1 that depends on the geometri- 
cal characteristics of the problem (Bruch, 1980). 
The w variable is numerically expressed with 
finite difference linear equations for the nodes 
of the mesh that discretize the flow zone under 
study. Applying this method, it is possible to 
determine the position of the upper flow line 
in a homogenous medium because it generates 
values of zero (w = 0) for the nodes with a total 
water pressure of zero (this means that these 
points are at atmospheric pressure). The form of 
the expression that describes w depends on the 
location of the nodes in the mesh. Consequently, 
the different expressions for obtaining w are 
described in the following paragraphs (Bruch, 
1980). For nodes inside the flow region, w is 
determined as follows: 


2,72 
A di - AA, | 
1,] 2 2 
2(A: + A5) 


1 m+l m 1 m+l m 
aer +, ¿+ A + w',x) = ] (22) 
dd y 


As well: 


o (0,01, + aut; e 0) (23) 


ij 


where a is the relaxation factor, max [...) is the 
maximum absolute value between the two data 
values of equation (23), and the remaining vari- 
ables are defined in the previous sections. For 
the nodes located in the upstream equipotential 
boundary, the variable w is defined as: 


1 2 
w= ¿a -L) (24) 


For the node located in the downstream 
equipotential boundary, w is calculated by: 


1 2 
w=>(L,-L) (25) 
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where L, is the height of the upstream water 
level (UWL), L, is the height of the downstream 
water level (DWL), and L is the height from a 
point of reference (for example, the foundation 
soil of the dam) to each one of the nodes accor- 
ding to the studied case. 

For nodes located in the lower flow bound- 
ary (base of the dam), the solution for w is: 


2 2 2 


where X, represents the length of the base of the 
dam in X direction, and x is the location of each 
node in X direction and 0 (origin) is the toe of 
the dam slope. 


Solutions to the Problems in Heterogeneous 
Media Using the ExtendedPressure Method 
(Brezis et al., 1978; Bardet € Tobita, 2002) 


In this method, the hydrostatic pressure p is 
assumed as the unknown quantity instead of 
the potential function q. Moreover, the Extended 
Pressure concept (Brezis et al., 1978; Bardet éz 
Tobita, 2002) is used, which modifies Darcy's 
law as follows: 


=-k-[grad(p)+ H.(p)grad(y)] (27) 


where e is the separation (Ay) between two 
nodes in the Y direction. Additionally: 


1  ifpze 
H (p)= 28 
¿p) a (28) 


The solution to equation (27) is determined 
similarly to obtaining the solution of Laplace's 
equation using the SOR method (equation 16) 
by introducing the relaxation factor a. Again, 
the expressions that govern the behavior at 
each node depend on their locations in the 
flow domain. Thus, the equations that define 
p are described below (Bardet €z Tobita, 2002). 
To evaluate the behavior of the nodes located 
inside the flow region, the hydrostatic pressure 
p is defined as: 


pr =(1-0)p", 


1 m+l m+1 
+ q Pija Piojo Pia Pia 


y 1 pul ea] 5) 
ñ 2 


For the nodes located in an impervious 


+A, 


lower flow boundary, the hydrostatic pressure 
is determined as follows: 


qe = (1- a ) cade m 


A, 
a apt jA + Pia j + e ) + on (30) 


In contrast to the Baiocchi's method, the 
value that is assigned to the nodes at the equi- 
potential boundaries of the mesh is equal to the 
hydrostatic pressure regarding the upstream 
or downstream water levels (UWL or DWL, 
respectively) for each specific case when using 
the Extended Pressure method. 

For the nodes located at the transition 
between two materials with different perme- 
ability's, the following expression is used to 
calculate the pressure: 


m+l 


po =0-o)p;, 


m+1 m 


1 m+ m 
as (2p Pi e +(1+ y Mp. 1 +P; 2 Pa 


PijtPija 
31 
m0 J 1) 


For a node located in the transition zone 


+A, 


y Po Pal 
2 


and in the impervious lower boundary, the ap- 
plicable equation is: 


pi ap! 


m m+l 


Pij+Pija 
2 


m+l m+l 


L Pra tU Pra + Pia 
1+y 


+ 


+AyH, 


Dd 


| (32) 
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where y = k,/k,. Here, k, is the permeability of 
the first material and k, is the permeability of 
the second material. 


Practical Applications 


The Initial Conditions using the Baiocchi's 
and Extended Pressure Methods 


One advantage of working with equations in 
the nodes of a mesh is that they can be pro- 
grammed without major complications using 
spreadsheets without more sophisticated coding 
or programming in other types of software. The 
equations of each node are independently and 
numerically solved by only considering the data 
of the adjacent nodes. In consequence, it is not 
necessary (as in the finite element method) to 
study the behavior of each element, and then 
with the assembly to solve complex linear equa- 


tion systems when evaluating global behavior. 
This process involves handling and storing a 
large number of matrices. In the following para- 
graphs, the use of spreadsheets is illustrated by 
using the equations of the SOR method to solve 
two free surface (unconfined flow) problems in 
dams. 

Figure 2 shows an unconfined flow problem 
in a homogenous and isotropic dam. In the same 
figure, the locations of some of the nodes and 
boundaries with the equations or conditions 
that govern their behavior are shown. These 
parameters are listed in Table 1 for the Baioc- 
chi's and Extended Pressure methods. 

Similarly, Figure 3 exemplifies a flow prob- 
lem in a dam that is constructed of materials 
with different permeability's. The boundary 
conditions and their governing equations for 
this case when applying the Extended Pressure 
method are provided in Table 2. 


Figure 2. Approach of a surface free problem (unconfined flow) through a homogenous and isotropic earth dam. 


Table 1. Baiocchi's and Extended Pressure methods, parameters. 


Boundary Baiocchi Extended Pressure 
AB Eq. (24) p=L,-L 
AE Eq. (26) Eq. (30) 
DE Eqs. (22) and (23) Eq. (29) 
FG Eq. (25) p=L,-L 
GH w=0 p=0 
J Egqs. (22) and (23) Eq. (29) 
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Figure 3. Approach of a surface free problem (unconfined flow) through a dam composed of materials with different 


permeability's. 


Table 2. Boundary conditions and their governing equations. 


Boundary Extended Pressure Boundary Extended Pressure 
AB p=L,-L GH p=0 
AE Eg. (30) I Eq. (29) 
DE Eq. (29) J Eg. (32) 
FG p=L,-L K Eg. (31) 


It is recommended that the nodes in the 
downstream slope of the dam or earth structure 
(boundaries FGH, indicated with open circles in 
Figure 2 and Figure 3, respectively) should be 
placed outside the flow zone (the ghost nodes) 
(Wang €: Anderson, 1982). This process allows 
the upper flow line to develop fully inside 
the flow zone rather than being forced to end 
exactly at the DWL (downstream water level). 
Otherwise, it is possible that the free discharge 
surface will not be well defined. 

In the problems discussed herein, a relaxation 
factor of a.=1.7 was assumed, which is the recom- 
mended value for significantly optimizing and 
reducing the calculation time (Salmasi € Aza- 
mathulla, 2013). Besides, the analyses are con- 
ducted in the following stages: 1) calculation of 
the upper flow line position, 2) evaluation of the 
potential function q for drawing up equipotential 
lines, 3) estimation of the stream function values 
Y for drawing up the flow lines, 4) calculation of 


the hydraulic gradients, and 5) evaluation of the 
seepage flux (as described below). 


Flow Problem Solution for a Homogenous 
Dam 


The flow of water through a homogenous 
earth dam that was built on impervious rock 
was studied. Figure 4 shows the geometrical 
characteristics and the mesh that were used for 
implementing the SOR method, with a node 
spacing or distance of Ax = Ay = 0.20 m. The 
studied dam has 2:1 slopes upstream and down- 
stream with assumed water levels of UWL= 9.0 
m (upstream) and DWL = 2.0 m (downstream). 

Along with a spreadsheet, Figure 5 presents 
the results that were obtained in stage 1) when 
deducing the position of the upper flow line. 
The spreadsheet was cut so to appreciate the 
results at the toe of the dam and near its crest (in 
part of the upstream slope). The values in each 
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Figure 4. Geometric characteristics and mesh used for implementing the SOR method in a homogenous earth dam that was 


built on impervious rock. 


of the cells represent the solutions that were 
obtained by means of the Baiocchi's (Figure 5a) 
and Extended Pressure (Figure 5b) equations, 
which are indicated in Section “The initial 
conditions using the Baiocchi's and Extended 
Pressure methods”. 

Once the upper flow line is defined, the 
boundary conditions are modified for continu- 
ing the calculation process. First, the potential 
function q is determined for drawing up the 
equipotential lines. Second, the values of the 
stream function Y are obtained for drawing up 
the flow lines. In this way, in stage 2), the respec- 
tive values of total pressure h (Equation (3)) are 
assigned in the nodes located at the boundaries 
of water infiltration (upstream), the discharge 
(downstream) area, and in the upper flow line to 
evaluate the potential function q. Whereas, the 
nodes of the lower flow boundary (impervious 
base of the dam) are governed by Equation (17). 
In stage 3) considering an analogy of Y with 
q, it is possible to use the same equations for 
calculating the stream function Y that were used 
for calculating the potential function q. Here, 
the boundary conditions are modified again 
according to the following situations: a) in the 
lower flow line (impervious boundary), a value 
of zero is assigned to all nodes (Y, = 0) and b) 
the value assigned to the upper flow line repre- 
sents the flux that passes between the upper and 
lower flow lines, which corresponds to the total 
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seepage flux through the dam (Y, = q,,,) (Bardet 
£z Tobita, 2002; Bosch éz Davis, 1997). Figure 6a 
exhibits the position of the upper flow line that 
was obtained using the Baiocchi's and Extended 
Pressure techniques (based on the SOR method). 
For comparison purposes, this figure also shows 
the position of the upper flow line that was cal- 
culated with the Kozeny”s parabola (a rigorous 
method proposed by Kozeny in 1931) and the 
finite element method (an approximate numeri- 
cal method that was applied here through the 
SEEP/W algorithm, Geo-Slope International 
Ltd., 2008). 

As revealed in Figure 6a, the upper flow 
lines calculated with the procedure proposed 
here using the SOR method were similar to 
those determined by the other methods. In 
particular, the Extended Pressure technique 
matches the solution that was obtained by the 
finite element method. Figure 6b provides the 
numerically calculated flow net based on the 
SOR method, which considers the position of 
the upper flow line that was obtained with the 
Extended Pressure method. Additionally, Figure 
6b confirms that the SOR method can be used to 
define the flow net accurately. 


Flow Problem Solution for a Dam Composed 
of Materials with Different Permeability 's 


Here, the water flow through a dam compo- 
sed of distinct materials was studied. Figure 
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a) Values of w computed with Baoicchi's solution 


17.33 | 17.01 | 16.69 
8.80 A 9.47 | 19:14 | 1880 | 18.47 | 1814 | 17.82 | 17.49 
9.00 Ro 19.66 | 19.32 | 18.98 | 18.64 
9.20 20.88 | 20.53 | 20.18 | 19.84 
9.40 23.21 | 2285 | 2249 | 2214 | 21.78 | 21.43 | 21.07 
9.60 23.44 | 23.07 | 22.71 | 2234 
9.80 24.78 | 2440 | 24.03 | 23.65 
10.00 26.46 | 25.77 | 25.39 | 25.01 
x(m) 0.00 18.80 | 19.20 | 1960 | 20.00 
Boundary values (Eq. 24 
== y (Eq-24) Internal values (Eq. 22 and 23) 
Boundary values (Eq. 26) 
b) Values of p calculated with Extended Pressure method 
B [e D E F G H I J K AU AV AW AX AY AZ BA BB 
¡) 
y (m) 
000 | ===-=-- 7 | 
0.20 DD e | 0.00 [| 0.00 
0.40 | 000 [ o.00 | 0.00 
0.60 > | 0.00 | 0.00 0.00 | 0.00 
0.80 | 0.00 | 0.00 0.00 0.00 | 0.00 
1.00 — e ——— [ 000 [| 0.01 | 0.01 0.00 0.00 | 0.00 
1.20 0.08 0.04 | 0.02 0.01 0.01 | 0.00 
1.40 040 | 029 0.16 0.09 0.05 0.03 0.02 0.01 
8.60 7.16 7.10 7.05 699 | 99 6.86 6.79 | 6.73 
8.80 Ej 34 | 728 | 722 | 716 | 710 | 703 | 696 
9.00 1 8.04 | 8.07 7A 3 [58 7.52 7.46 7.40 7.33 vay 7.20 
9.20 8.24 827 | 8.31 7.87 7.82 7.76 7.70 | 7.64 7.58 7.51 7.44 
9.40 8.44 | 8.48 8.52 | 8.55 8.11 8.06 8.00 7.94 | 7.88 7.82 | 7.75 7.69 
9.60 8.64 8.68 | 8.72 8.76 | 8.79 8.36 8.30 8.25 8.19 | 8.13 8.06 8.00 7.93 
9.80 8.85 | 8.89 8.93 8.97 9.01 9.04 8.61 8.55 8.49 8.43 | 8.37 8.31 8.24 8.18 
10.00 900 | 905 | 909 | 914 9.18 | 9.22 9.25 | 9.29 8.85 | 8.80 8.74 | 8.68 | 8.62 8.56 | 849 | 8.43 
x(m) | 0.00 0.40 | 0.80 | 120 160 | 200 | 240 | 280 | 1720 | 17.60 | 18.00 | 18.40 | 18.80 | 1920 [1960 | 20.00 
Boundary values (hydrostatic pressure 
[1] y (hy Pp ) Internal values (Eq. 29) 


Boundary values (Eq. 30) 


Figure 5. Values of w and p that were calculated at the nodes of boundaries and inside the flow zone (upstream slope) by 
applying; a) the Baiocchi's solution, and b) the Extended Pressure method. 


7 shows the geometric characteristics and the UWL = 9.0 m (1.0 m below the crest of the dam) 
mesh that was used for implementing the SOR and upstream and downstream slopes of 2:1 
method with a node spacing or distance of for the dam. The core of the dam is composed 
Ax =0.20 and Ay =0.10 m, respectively. Likewise, of a low permeability material (k,) relative to 
this figure shows an upstream water level of the adjacent transition materials (k, and k,). The 
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a) 


UWL 
Y 


- =--- Kozeny's solution 
—— Finite element 

— Extended Pressure method 
== Baiocchi's solution 


___ Upper flow line 

(Extended Pressure method) 
—- Equipotential lines 
—— Flow lines 


Figure 6. Study of the flow of water through a homogenous and isotropic earth dam: a) comparison of the upper flow lines 


obtained by using different methods, and b) the flow net numerically calculated with the Extended Pressure method (based on 


ratio between the permeabilities of the core 
material and the upstream transition material is 
k,/ k,=0.1, and the ratio between the core mate- 
rial permeabilities and the downstream transi- 
tion material is k,/ k,=10. At the toe of this dam, 
a horizontal sand filter with high permeability 
and a length of 10.0 m was placed, which relie- 
ves water pressure to avoid or mitigate erosion 
problems at this part of the dam. 

This practical application is solved similarly 
to the previous application as follows: stage 1) 
deduction of the upper flow line position with 
the Extended Pressure method, stage 2) evalua- 
tion of the potential function q for drawing the 
equipotentials lines, and stage 3) calculation of 
the stream function values Y for drawing the 
flow lines. Moreover, once the water pressure 
variations in the dam are obtained, in order 
to complete the flow problem solution the 
hydraulic gradients are determined in stage 4) 
and the seepage flux is calculated in stage 5) 
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as indicated below. To evaluate the hydraulic 
gradient (resulting magnitude) inside the flow 
zone, the following expression can be applied in 
terms of finite differences (Budhu, 2000): 


Ej] [Ea 


¡(res),, = e + na 


2 


(33) 


where ¡(res) is the hydraulic gradient value 
(resulting magnitude) in each node of the mesh. 

The flow rate calculation is determined by 
drawing a vertical line or plane across the flow 
region and using the pairs of nodes over this 
line that are located in transverse direction to 
the flow. Thus, the flow rate is obtained from 
the following equation (Budhu, 2000): 


KA (E 
Uos = => d. AT b., ¡ (S4) 
2A yl 15 ») 


ha 
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Figure 7. Geometric characteristics and mesh used for implementing the SOR method in a dam composed of different materials 


and a high permeability sand filter at the toe of the downstream slope. 


where: 
1 = total flow rate across the flow zone. 
k. = permeability in the X direction. 


D = nodes located in the lower line. 

E = nodes located in the upper line. 

Ax = distance between two nodes in the X 
direction. 

Ay = distance between two nodes in the Y 
direction. 


Figure 8 shows the spreadsheet that is used 
to solve the problem. The values of the cells cor- 
respond with the hydrostatic pressure (p) that 
was calculated for the nodes of the impervious 
boundary, the filter and inside the flow zone 
when applying the Extended Pressure method 
(the spreadsheet is cut off to show a portion of 
the results). 

Figure 9a shows the flow net that was nu- 
merically obtained with the SOR method in 
the analyzed dam. In this case, the upper flow 
line was calculated with the Extended Pressure 
technique. For comparison purposes, Figure 9b 
presents the upper flow line that was obtained 
using the finite element method through the 
SEEP/W algorithm (Geo-Slope International 
Ltd., 2008). When comparing the results ob- 
tained with SOR and FEM, similarities were 


observed. In the zone of the material with lower 
permeability, a greater loss of hydraulic head 
occurs. However, due to the law of continuity 
for the steady-state flow, the total flow rate that 
traverses the three flow zones is the same q,., 
= 0.65 m*/s. This flow rate was calculated in 
terms of finite differences with Equation (34). 
Similarly, the classical theory notes (Cedergren, 
1989) cases where the flow domain consists of 
two or more portions with different perme- 
abilities (each one consists of a homogenous and 
isotropic soil), in which the flow net is distorted 
at the boundaries between contiguous materials 
(Figures 9a and 9b) so that the same amount of 
flux passes through both sides of the boundary 
between two flow lines. Thus, Figures 9a and 
9b point out that the permeability is lower for 
the wider flow channels and greater for the nar- 
rower flow channels (Flores, 2000). 

Finally, Figure 10 provides the hydraulic 
gradient (resulting magnitude) in the flow zone 
that was calculated using finite differences with 
Equation (33). 


Conclusions 
Here, an alternative for analyzing free surface 


problems was proposed. This alternative was 
established on the numerical solution of finite 


ISSN 0187-8336 + 


Tecnología y Ciencias del Agua, vol. VI, núm. 5, septiembre-octubre de 2015, pp. 43-58 


López-Acosta £ González-Acosta, Study of Water Flow in Dams using Successive Over-relaxation 


1 171 172 

0.01 [ ooo [ ooo | o.w00 [ ooo [ ooo | o00 [ooo [ o00 ]ow00 [ooo [ooo | o00 [ ooo [ o.00 | o.00 
y (m) 0.01 | oo | ooo | o00 | ow00 | o00 | ow00 [ooo [ o00 |o00 | ow00 | o00o | o00 | oo00 | o00 | 0.00 
82 0.02 | 001 | o01 | o00 | ow00 | o00 | oo00 | o00 [ooo | oo | o00 | oo | oo00 | oo00 | 000 | 000 
83 0.04 | oo | oo [ oo | oo | ooo | o00 | ooo | o00 [o00 | oo [ooo [ o00 | o00 [000 [ 0.00 
8.4 0.08 | 0.05 [| o03 | oo0z | o01 | oo | ooo [ooo | ooo [ooo | o00 | ooo [ooo | o00 [ o00 [0.00 
8.6 0.19 | 014 | o09 | o05 | o03 | oo | oo1 | oo01 | o00 [ooo | o00 | o00 | o00 | o00 | o00 | 0.00 
87 025 | 020 [| o14 [ 009 | o05 | oo | o02 | oo | oo [o00 | o00 | o00 [o00 | o00 [ o00 [ 0.00 
8.8 031 | 025 [| o19 [ o14 | oo0s | oos | o03 | oo | oo1 [oo | o00 | ooo [ooo | o00 [| o00 [ 0.00 
89 036 | o3o | o24 | o18 | o13 | o0s [ o04 | o03 [ oo [ o01 | oo | o00 | o00 | o00 | ow00 | 000 
9 o42 | 035 | o29 | o23 | o17 | o11 | o07 | o0% | oo | oo01 | oo | o00 | o00 | oo00 | o00 | 000 
9.1 0.47 | oso [| o34 | o27 | o21 | o15 [| o10 | o06 | oo3 [o02 | oo01 | oo | o00 | o00 [ o00 | 000 
92 052 | 045 | o38 | o31 | o25 | o18 [| o13 | oos | oos [o03 | oo | oo01 | oo [ oo00 | o00 | 000 
93 os5s | oso | o“ | o35 | o28 | o21 | o15 | o.10 | o06 [ o04 | oo | oo01 | o01 | o01 | 000 | 000 
9.4 o.63 | 054 | 046 | o38 | o30 | o23 | o.17 | o.12 | oos [ o0s | oo03 | oo | o01 | oo1 | 000 | 000 
95 o.68 | 059 | 09 | oso | os32 | o25 | o19 | o.13 | oo09 [ o06 | o0% | oo | o01 | oo1 | 001 | 000 
96 0.73 | 063 | 053 | 04 | 033 | o25 | 019 | o14 | o1o [006 | oo | oo3 | o02 | oo | oo [ 000 
97 079 | 067 | 055 | 043 | 033 | o24 | o18 | o13 | o10 [o07 | oo0s | oo [ 002 | o01 [| o01 | 001 
9.8 os | 072 | o58 | 043 | o3o | o21 | o15 | 0.12 | o09 [ o07 | oos | oo03 | o02 | o01 [ o01 | 001 
99 093 | 078 | 062 [ 043 | 022 | o14 | o1m1 | o09 | o07 [| o05 | o0% | oo3 [ 002 | oo [ oo1 | 000 

10 1.02 0.87 0.68 0.43 
x(m) | 332 | 334 | 356 | 338 | 34 | 342 | 344 | 346 | 348 | 35 352 | 354 | 356 | 358 | 36 36.2 

E] Boundary values (filter hydrostatic pressure) Intemál vábies (Eg. 29) 


Boundary values (Eq. 30) 


Figure 8. Values of the hydrostatic pressure (p) that were calculated for the nodes of the impervious boundary, in the filter and 


inside the flow zone by applying the Extended Pressure method. 


a) Aro = 0.65 m*/s 
A 
| ___Upper flow line 
gut (Extended Pressure method) 
AAA A te —— Equipotential lines 
—— Flow lines 


SN PO 
rs TT 
b) 
Upper flow line 


(finite element) 
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, núm. 


CA RIA 


7, 


Figure 9. a) Flow net that was calculated with the SOR method in a dam composed of different materials and a horizontal filter 
at the toe of the downstream slope, and b) the upper flow line obtained by the finite element method. 
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UWL 
Ze 


___ Upper flow line 
(Extended Pressure method) 


—— Resulting hydraulic gradient 


Figure 10. Resulting hydraulic gradient value that was calculated by the SOR method. 


difference equations using the SOR (Successive 
Overrelaxation) method. The implementation of 
two techniques based on the SOR method were 
shown, including the Baiocchi's solution and the 
Extended Pressure method, which utilize the 
Gauss-Seidel iterative process. The equations 
of the SOR method can be easily enabled using 
spreadsheets for the solutions of different types 
of flow problems of variable complexity (such 
as free surface problems) without requiring so- 
phisticated programing or specialized software. 
Using the SOR method, the equations of each 
node are numerically solved independently and 
automatically by considering the data of the 
adjacent nodes. Thus, it is not necessary (as in 
other numerical methods such as the FEM meth- 
od) to study the behavior of each element and 
afterward with the assembly to solve complex 
systems of linear equations in order to evaluate 
the global behavior, which involves the storage 
and handling of a large number of matrices that 
increase the calculation time. The applicability 
of the method was exposed by the unconfined 
water flow analysis through homogenous and 
heterogeneous earth dams. By using the SOR 
technique, the variations of the potential func- 
tion q and the stream function Y were calculated 
for numerically drawing the flow net. In addi- 
tion, the resulting magnitude of the hydraulic 
gradient, the total flow rate and the position 
of the upper flow line were estimated with the 
Baiocchi's and Extended Pressure techniques. 


From the analyses performed, the Extended 
Pressure method was shown to provide results 
with better approximation than the Baiocchi's 
solution. In particular, the Extended Pressure 
method provides practically equal results to 
those calculated with FEM. One additional 
advantage that was observed in the Extended 
Pressure method is that it does not require a 
very refined discretization of the flow domain. 
In consequence, this method can be used to 
obtain good results even when the separation 
between the nodes of the mesh is large (a coarse 
grid). The obtained results demonstrate that 
the techniques implemented in this paper are 
simple and easily applied. Furthermore, these 
techniques provided similar results to those 
obtained with the currently popular numerical 
methods, such as the finite element method. 
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